#include "Mesh.h"

#define ACCURACY 1e-12
#define MIN(a, b) ((a) > (b) ? (b) : (a))
#define MAX(a, b) ((a) < (b) ? (b) : (a))
#define SIGN(x) ((x) > 0 ? 1 : -1)

using namespace CG;
using namespace LA;

// 创建顶点，并返回指针
vertex3f *Mesh::CreateVertex(point3f p, bool testSame)
{
	// 如果开启检测，会禁止构造重合的顶点
	if (testSame)
	{
		for (auto it = m_vList.begin(); it != m_vList.end(); ++it)
		{
			// 返回重合顶点
			if (length((*it)->pos - p) < ACCURACY)
				return *it;
		}
	}
	vertex3f *vp = new vertex3f({p, nullptr});
	m_vList.push_back(vp);
	return vp;
}

// 创建三角面，可选自动衔接
face *Mesh::CreateFace(vertex3f *v0, vertex3f *v1, vertex3f *v2, bool autoTwin)
{
	// 创建三角形面
	vertex3f *vList[] = {v0, v1, v2};
	face *f = CreateFace(vList, 3, autoTwin);

	// 如果之前是空网格，现在设为三角形网格
	if (m_type == MeshType::empty)
		m_type = MeshType::triangular;
	return f;
}

// 创建四边形面，默认自动判定衔接相邻面
face *Mesh::CreateFace(vertex3f *v0, vertex3f *v1, vertex3f *v2, vertex3f *v3, bool autoTwin)
{
	// 判断四点是否共面，如果不共面，转为建立两个三角面
	if (!onFace(v0->pos, v1->pos, v2->pos, v3->pos))
	{
		// 注意创建顺序：先创建左下角的三角形，然后创建右上角的三角形
		// 这个顺序的目的是当从左向右，从下向上创建面时，不会出现流形情况
		CreateFace(v0, v1, v3);
		CreateFace(v1, v2, v3);
		return nullptr;
	}

	// 创建四边形面
	vertex3f *vList[] = {v0, v1, v2, v3};
	face *f = CreateFace(vList, 4, autoTwin);

	// 如果之前是空网格，现在设为四边形网格
	if (m_type == MeshType::empty)
		m_type = MeshType::rectangular;
	// 如果之前是三角网格，现在就是混合网格
	else if (m_type == MeshType::triangular)
		m_type = MeshType::mixed;
	return f;
}

// 创建多边形面
face *Mesh::CreateFace(vertex3f **vList, int n, bool autoTwin)
{
	// 有一个顶点为空，就退出
	for (int i = 0; i < n; i++)
		if (vList[i] == nullptr)
			return nullptr;

	face *f = new face;
	half_edge **eList = new half_edge *[n];

	// 创建半边
	for (int i = 0; i < n; i++)
		eList[i] = new half_edge;

	// 初始化半边
	for (int i = 0; i < n; i++)
	{
		half_edge *e = eList[i];

		e->v = vList[i];
		e->incFace = f;
		e->next = eList[(i + 1) % n];
		e->prev = eList[(i + n - 1) % n];
		e->twin = nullptr;
	}

	// 面持有任意一个半边
	f->e = eList[0];

	// 默认自动衔接面
	if (autoTwin)
		CreateTwin(f);

	// 如果顶点没有半边，就赋予半边，放在这里，因为如果顶点没有半边，不需要计算邻面
	for (int i = 0; i < n; i++)
	{
		if (vList[i]->e == nullptr)
			vList[i]->e = eList[i];

		// 推入半边
		m_eList.push_back(eList[i]);
	}
	delete[] eList;

	// 最后推入面
	m_fList.push_back(f);

	return f;
}

// 为一个面创建连接，会选择面的顶点相邻的面尝试建立连接
void Mesh::CreateTwin(face *f)
{
	std::deque<vertex3f *> vList = VertexOnFace(f);

	// 遍历整个面表，逐一尝试建立 twin 连接
	// 优化：不再遍历面表，而是选择与顶点相邻的面尝试建立 twin 连接；当面的数量较大时具有优势
	// 再优化：由于三角形顶点两两共用面，只需要计算两个顶点的邻面（有问题，不能只计算两个顶点，需要三个顶点都计算）
	// 注意要先取出所有邻面再建立连接，防止出现与自己连接的情况
	std::deque<face *> fList0 = AdjacentFace(vList[0]);
	std::deque<face *> fList1 = AdjacentFace(vList[1]);
	std::deque<face *> fList2 = AdjacentFace(vList[2]);
	for (auto it = fList0.begin(); it != fList0.end(); ++it)
		StitchFace(f, *it);
	for (auto it = fList1.begin(); it != fList1.end(); ++it)
		StitchFace(f, *it);
	for (auto it = fList2.begin(); it != fList2.end(); ++it)
		StitchFace(f, *it);
}

// 删除顶点
bool Mesh::DeleteVertex(vertex3f *v)
{
	if (v == nullptr)
		return false;

	// 获得所有相邻的面，将它们删除
	std::deque<face *> fList = AdjacentFace(v);
	for (auto it = fList.begin(); it != fList.end(); ++it)
		DeleteFace(*it);

	// 最后将 v 孤立
	v->e = nullptr;

	return true;
}

// 删除半边和相邻的面
bool Mesh::DeleteHalfEdge(half_edge *e)
{
	// 要求网格是三角形网格，要求对边必须存在
	if (m_type != MeshType::triangular || e == nullptr || e->twin == nullptr)
		return false;

	// 将 v2 合并到 v1
	vertex3f *v1 = e->v;
	vertex3f *v2 = Destination(e);

	// 获得需要操作的边，e11 与 e21 要配对，e12 与 e22 要配对
	half_edge *e11 = e->prev->twin;
	half_edge *e12 = e->twin->next->twin;
	half_edge *e21 = e->next->twin;
	half_edge *e22 = e->twin->prev->twin;

	// 这里记录 v1 将要保存的半边，因为删除面的操作可能使得 v1,v2 孤立
	half_edge *tmpe = nullptr;
	// 检查特殊情况：即 e11-e12 或 e21-e22 配对，这意味着两边重合，将会被删除，需要重新指定配对方式
	if (e11->twin != e12)
		tmpe = e11;
	// 如果 e11-e12，那么 e12->next 将会被删除，此半边无效
	if (e12->twin != e11)
		tmpe = e12->next;
	if (e21->twin != e22)
		tmpe = e21->next;
	if (e22->twin != e21)
		tmpe = e22;

	// 在删除面之前要先检查配对方式。正常的配对方式是 e11-e21, e12-e22
	// 但是，如果 e11,e12 或 e21,e22 原本就配对，意味着 e11,e12 或 e21,e22 将会被删除
	// 此时需要改变配对方式为 e21-e22 或 e11-e12
	int twinType = 0;
	if (e11->twin == e12)
		twinType = 1;
	if (e21->twin == e22)
		twinType = 2;
	// 有可能出现 e11-e12, e21-e22 的情况，但是这种情况在前面会造成 tmpe == nullptr，从而退出，这里不需要再考虑

	// 记录目标是 v1 的边数
	int count = 0;
	std::deque<half_edge *> eList = AdjacentHalfEdge(v2);
	for (int i = 0; i < eList.size(); i++)
	{
		if (Destination(eList[i]) == v1)
			count++;

		// 如果目标是 v1 的边数多于 1 个，说明在 v1-v2 构成的边上，存在多于 2 个半边
		// 但我们这里只删除 2 个半边，因此剩余的半边会叠合起来
		// 并且我们删除了一个顶点，之后对剩余半边进行收缩时就不会删除顶点，这里立即返回
		DEBUG_TYPE(count < 2, DEBUG_ERROR_EDGE_MESH, false)
	}

	// 以 v2 为起点的半边，将它们的起点都改为 v1
	for (int i = 0; i < eList.size(); i++)
		eList[i]->v = v1;

	// 令 v1 孤立，这样删除面时不需要再重新赋予 v1 半边
	v1->e = nullptr;

	// 进行删除操作，对于边收缩而言，会删除 1 到 2 个面
	// 第一个面是从 v1 指向 v2 的半边对应的面，这个面一定存在，因为 e 存在；第二个是该半边的反向边对应的面，不一定存在
	// 注意先删除 twin，不然 e 先被删除就成为 nullptr 了
	DeleteFace(e->twin->incFace);
	DeleteFace(e->incFace);

	// 令 v2 孤立，同时赋予 v1 新的半边
	v2->e = nullptr;
	v1->e = tmpe;

	// 基本配对方式
	if (twinType == 0)
	{
		e11->twin = e21;
		e21->twin = e11;
		e12->twin = e22;
		e22->twin = e12;
	}
	// e11-e12 的情况，配对 e21-e22
	else if (twinType == 1)
	{
		e21->twin = e22;
		e22->twin = e21;
	}
	// e21-e22 的情况，配对 e11-e12
	else
	{
		e11->twin = e12;
		e12->twin = e11;
	}

	return true;
}

// 删除面
bool Mesh::DeleteFace(face *f)
{
	if (f == nullptr)
		return false;

	// 获得三条半边，断言它们都非空
	std::deque<half_edge *> e = HalfEdgeOnFace(f);
	for (int i = 0; i < e.size(); ++i)
		DEBUG_ASSERT(e[i] != nullptr, DEBUG_NULL_HALFEDGE, false)

	// 调整每条半边起点保存的半边
	for (int i = 0; i < e.size(); i++)
	{
		// 获得半边的起源点，和以 v 为起点的半边
		vertex3f *v = e[i]->v;
		DEBUG_ASSERT(v != nullptr, DEBUG_NULL_VERTEX, false)

		// 当起点保存的半边是此半边时，就修改保存的半边。因此可以提前修改，这里就不会修改
		if (v->e == e[i])
		{
			// 删除面时，如果面上的顶点已经孤立，说明出现了非流形网格
			// 此时邻边的数量为零，应当重点注意
			std::deque<half_edge *> eList = AdjacentHalfEdge(v);

			// 如果 v 不止一条半边，就赋予新的半边，否则 v 将会孤立
			for (auto it = eList.begin(); it != eList.end(); ++it)
			{
				// 赋予不同的半边
				if ((*it) != e[i])
				{
					v->e = *it;
					break;
				}
			}
		}
	}

	// 删除每条半边
	for (int i = 0; i < e.size(); i++)
	{
		// 首先将每条半边的反向边的反向边设为空
		half_edge *te = e[i]->twin;
		if (te != nullptr)
			te->twin = nullptr;

		// 如果此半边的起点仍然持有此半边，那么起点将会孤立，将该点保有的半边设为空
		if (e[i]->v->e == e[i])
			e[i]->v->e = nullptr;

		// 将半边保有的顶点设为空
		e[i]->v = nullptr;
	}

	// 将面保有的半边设为空
	f->e = nullptr;

	return true;
}

// 清除孤立点
bool Mesh::DeleteIsolated(CGElement type)
{
	bool ifDelete = false;

	switch (type)
	{
	case CGElement::vertex:
	{
		// 找到所有孤立点删除
		for (auto it = m_vList.begin(); it != m_vList.end();)
		{
			// 如果不存在以它为起点的半边，说明此顶点孤立，直接删除
			if ((*it)->e == nullptr)
			{
				it = m_vList.erase(it);
				ifDelete = true;
			}
			else
				++it;
		}
		break;
	}
	case CGElement::half_edge:
	{
		// 找到所有孤立边删除
		for (auto it = m_eList.begin(); it != m_eList.end();)
		{
			// 如果不存在半边起点，说明孤立
			if ((*it)->v == nullptr)
			{
				it = m_eList.erase(it);
				ifDelete = true;
			}
			else
				++it;
		}
		break;
	}
	case CGElement::face:
	{
		// 找到所有孤立面删除
		for (auto it = m_fList.begin(); it != m_fList.end();)
		{
			// 如果不存在半边，说明孤立
			if ((*it)->e == nullptr)
			{
				it = m_fList.erase(it);
				ifDelete = true;
			}
			else
				++it;
		}
		break;
	}
	}

	return ifDelete;
}

// 获得网格的包围盒
bounding_box Mesh::GetBox()
{
	m_box.left = 0;
	m_box.right = 0;
	m_box.top = 0;
	m_box.bottom = 0;
	m_box.zFar = 0;
	m_box.zNear = 0;

	// 空网格返回空
	if (VertexNum() == 0)
		return m_box;

	// 取一个点设置
	point3f pos = m_vList.front()->pos;
	m_box.left = pos.x;
	m_box.right = pos.x;
	m_box.top = pos.y;
	m_box.bottom = pos.y;
	m_box.zFar = pos.z;
	m_box.zNear = pos.z;

	// 遍历网格点计算
	for (int i = 0; i < VertexNum(); i++)
	{
		pos = m_vList[i]->pos;
		m_box.left = MIN(m_box.left, pos.x);
		m_box.right = MAX(m_box.right, pos.x);
		m_box.bottom = MIN(m_box.bottom, pos.y);
		m_box.top = MAX(m_box.top, pos.y);
		m_box.zFar = MIN(m_box.zFar, pos.z);
		m_box.zNear = MAX(m_box.zNear, pos.z);
	}
	return m_box;
}

// 是否有法向
bool Mesh::HasVertexNormal()
{
	return !m_nList.empty();
}

// 判断网格是否是凸包
bool Mesh::IsConvexHull()
{
	// 对每个顶点，遍历每一个面
	for (auto vit = m_vList.begin(); vit != m_vList.end(); ++vit)
	{
		// 要排除离散点
		if ((*vit)->e == nullptr)
			continue;

		// 判断顶点是不是在所有面的上方，如果在某个面的上面，说明不是凸包
		for (auto fit = m_fList.begin(); fit != m_fList.end(); ++fit)
			if (AboveFace((*vit)->pos, *fit))
				return false;
	}

	return true;
}

// 判断点是否在网格的凸包中
bool Mesh::IsInMesh(point3f p)
{
	// 遍历每一个面，p 如果在某个面的上面，说明 p 在网格外
	for (auto it = m_fList.begin(); it != m_fList.end(); ++it)
		if (AboveFace(p, *it))
			return false;
	return true;
}

// “缝合”两个面，返回是否出错
bool Mesh::StitchFace(face *f1, face *f2)
{
	// 先取得公共半边
	std::deque<half_edge *> e = CommonOnFace(f1, f2);

	// 如果有公共边就连接，一种情况是如果已经建立，但是不是两个面的公共边，直接报错
	if (e[0] != nullptr && e[1] != nullptr)
	{
		DEBUG_ASSERT(e[0]->twin == nullptr || e[0]->twin == e[1], DEBUG_EXIST_TWIN, false)
		DEBUG_ASSERT(e[1]->twin == nullptr || e[1]->twin == e[0], DEBUG_EXIST_TWIN, false)

		// 然后互相为 twin
		e[0]->twin = e[1];
		e[1]->twin = e[0];
		return true;
	}
	// 没有公共边也报错
	return false;
}

// 判断一个点是否在面上
bool Mesh::OnFace(point3f p, face *f)
{
	// 排除面是空的情况
	if (f == nullptr)
		return false;

	// 找到面的一个顶点
	point3f p0 = f->e->v->pos;
	vector3f v = p - p0;

	// 如果与法向正交，说明共面
	return (fabs(dot(v, FaceNormal(f))) < ACCURACY) ? true : false;
}

// 判断一个点是否在面的上方
bool Mesh::AboveFace(point3f p, face *f)
{
	// 排除面是空的情况
	if (f == nullptr)
		return false;

	// 取一个顶点计算夹角
	point3f p0 = f->e->v->pos;
	vector3f v = p - p0;

	// 这里精度需要仔细考虑，不然判断会有问题
	return (dot(v, FaceNormal(f)) > ACCURACY) ? true : false;
}

// Gouraud 插值：计算顶点的平均法向
vector3f Mesh::Gouraud(vertex3f *v)
{
	// 先获取所有相邻的面
	std::deque<face *> fList = AdjacentFace(v);

	// 累加然后规范化
	vector3f n = {0, 0, 0};
	for (auto it = fList.begin(); it != fList.end(); ++it)
		n = n + FaceNormal(*it);
	return normalize(n);
}

// 计算单位面法向
vector3f Mesh::FaceNormal(face *f)
{
	// 找到面的三个顶点计算法向
	point3f v0 = f->e->v->pos;
	point3f v1 = f->e->next->v->pos;
	point3f v2 = f->e->next->next->v->pos;
	return normalize(cross(v1 - v0, v2 - v1));
}

// 获得边的目标点
vertex3f *Mesh::Destination(half_edge *e)
{
	if (e == nullptr)
		return nullptr;
	e = e->next;
	return e->v;
}

// 获得一个面上的所有顶点
std::deque<vertex3f *> Mesh::VertexOnFace(face *f)
{
	half_edge *e = f->e;
	std::deque<vertex3f *> vList = {e->v};
	e = e->next;
	while (e != f->e)
	{
		vList.push_back(e->v);
		e = e->next;
	}
	return vList;
}

// 获得一个面上的所有半边
std::deque<half_edge *> Mesh::HalfEdgeOnFace(face *f)
{
	half_edge *e = f->e;
	std::deque<half_edge *> eList = {e};
	e = e->next;
	while (e != f->e)
	{
		eList.push_back(e);
		e = e->next;
	}
	return eList;
}

// 获取两个面的公共半边
std::deque<half_edge *> Mesh::CommonOnFace(face *f0, face *f1)
{
	std::deque<half_edge *> eList(2, nullptr);
	std::deque<half_edge *> eList0 = HalfEdgeOnFace(f0);
	std::deque<half_edge *> eList1 = HalfEdgeOnFace(f1);

	// 遍历两个面上的所有边
	for (int i = 0; i < eList0.size(); i++)
	{
		for (int j = 0; j < eList1.size(); j++)
		{
			half_edge *e0 = eList0[i];
			half_edge *e1 = eList1[j];

			// e0 的起点是 e1 的终点，e1 的起点是 e0 的终点时，说明是公共边，储存后直接返回
			if (e0->v == Destination(e1) && e1->v == Destination(e0))
			{
				eList[0] = e0;
				eList[1] = e1;
				return eList;
			}
		}
	}
	return eList;
}

// 获得相邻的所有顶点
std::deque<vertex3f *> Mesh::AdjacentVertex(vertex3f *v)
{
	std::deque<vertex3f *> vList;

	// 先获得以 v 为起点的半边，其目标点就是相邻顶点
	std::deque<half_edge *> eList = AdjacentHalfEdge(v);
	for (int i = 0; i < eList.size(); i++)
		vList.push_back(Destination(eList[i]));

	// 当半边非空时才考虑
	if (!eList.empty())
	{
		// 注意考虑 twin 为空的情况，此时少算一个顶点
		// 第一种情况，只按照顺时针方向获得了半边，此时第一个半边的前一个半边的起点是少算的顶点
		if (eList.front()->prev->twin == nullptr)
			vList.push_back(eList.front()->prev->v);
		// 第二种情况，按照逆时针方向获得了半边，此时最后一个半边的前一个半边的起点是少算的顶点
		else if (eList.back()->prev->twin == nullptr)
			vList.push_back(eList.back()->prev->v);
	}

	return vList;
}

// 获得以 v 为起点的半边
std::deque<half_edge *> Mesh::AdjacentHalfEdge(vertex3f *v)
{
	std::deque<half_edge *> eList;
	half_edge *e = v->e;

	// 先顺时针寻找
	while (e != nullptr)
	{
		eList.push_back(e);
		e = e->twin;

		// 反向边，然后找到下一条边
		if (e != nullptr)
			e = e->next;

		// 如果回到初始边，就退出
		if (e == v->e)
			break;
	}

	// 如果没有回到初始边，就逆时针寻找，排除孤立顶点
	if (e == nullptr && v->e != nullptr)
	{
		e = v->e->prev->twin;
		while (e != nullptr)
		{
			eList.push_back(e);
			e = e->prev->twin;
		}
		// 这里不需要判断回到初始边，因为如果一个方向不能返回，那么反过来也不能返回
	}
	return eList;
}

// 获得相邻的所有面
std::deque<face *> Mesh::AdjacentFace(vertex3f *v)
{
	// 先获得以 v 为起点的半边
	std::deque<half_edge *> eList = AdjacentHalfEdge(v);

	// 每条半边左边的面就是相邻的面
	std::deque<face *> fList;
	for (auto it = eList.begin(); it != eList.end(); ++it)
		fList.push_back((*it)->incFace);

	return fList;
}

// 获得相邻的所有面
std::deque<face *> Mesh::AdjacentFace(face *f)
{
	std::deque<face *> fList;
	std::deque<half_edge *> eList = HalfEdgeOnFace(f);
	for (int i = 0; i < eList.size(); i++)
		if (eList[i]->twin != nullptr)
			fList.push_back(eList[i]->twin->incFace);
	return fList;
}

// 取网格点沿着 pv 方向的投影多边形顶点
std::deque<vertex3f *> Mesh::MeshProject(vector3f pv)
{
	// 标记被照射的面，只需要保存一个面的指针，因为所有标记面都相邻
	std::deque<face *> fList;
	for (auto fit = m_fList.begin(); fit != m_fList.end(); ++fit)
	{
		// 如果法向与 pv 同向，说明被照射/投影正面
		vector3f n = FaceNormal(*fit);
		if (dot(pv, n) > 0)
		{
			(*fit)->identity = 1;
			fList.push_back(*fit);
		}
		else
			(*fit)->identity = 0;
	}

	// 检查被照射的面，标记相邻的未照射面
	std::deque<half_edge *> eList;
	for (auto fit = fList.begin(); fit != fList.end(); fit++)
	{
		// 检查面上的半边，记录边界边
		std::deque<half_edge *> eList0 = HalfEdgeOnFace(*fit);
		for (int i = 0; i < eList0.size(); i++)
			if (eList0[i]->twin->incFace->identity != 1)
				eList.push_back(eList0[i]);
	}

	// 遍历所有反向边，按照连接顺序推入顶点
	auto eit = eList.begin();
	std::deque<vertex3f *> vList;
	while (vList.size() < eList.size())
	{
		vList.push_back((*eit)->v);
		// 找到目标点对应的边
		for (auto it = eList.begin(); it != eList.end(); ++it)
		{
			if ((*it)->v == Destination((*eit)))
			{
				eit = it;
				break;
			}
		}
	}

	return vList;
}

// 将网格三角化，返回三角化后的面列表
std::deque<face *> Mesh::TriangulateFace(face *f)
{
	// 获得所有顶点，如果恰好有三个，直接返回
	std::deque<vertex3f *> vList = VertexOnFace(f);
	if (vList.size() == 3)
		return {f};

	// 否则先删除面
	DeleteFace(f);

	// 遍历顶点创建
	std::deque<face *> fList;
	for (int i = 1; i < vList.size() - 1; i++)
		fList.push_back(CreateFace(vList[0], vList[i], vList[i + 1]));
	return fList;
}

// 获得顶点的二次误差矩阵
Matrix Mesh::QEMatrix(vertex3f *v)
{
	Matrix Q(4, 4);

	// 获得顶点的相邻面，计算每个面的法向
	std::deque<face *> fList = AdjacentFace(v);
	for (int j = 0; j < fList.size(); ++j)
	{
		// 要计算平面方程 ax+by+cz+d=0
		vector3f normal = FaceNormal(fList[j]);
		double d = -dot(normal, {v->pos.x, v->pos.y, v->pos.z});
		Vector p = {normal.x, normal.y, normal.z, d};
		Q = Q + p * p.transpose();
	}
	return Q;
}

// 获得边上顶点收缩的最优解
Vector Mesh::MinVector(half_edge *e)
{
	// 线性求解器
	LinearSolver solver;

	// 获得两个端点
	vertex3f *v1 = e->v;
	vertex3f *v2 = Destination(e);

	// 求解极小化问题的解 v，通过顶点标记获得对应的矩阵
	Matrix Q = m_QList[v1->identity] + m_QList[v2->identity];

	// 获得线性系统的矩阵并求解
	Matrix Q1 = Q;
	Q1(3, 0) = Q1(3, 1) = Q1(3, 2) = 0;
	Q1(3, 3) = 1;
	Vector b = {0, 0, 0, 1};
	Vector V = solver.solve(Q1, b, true);

	// 获得记录状态，如果状态错误，说明求解有问题，改为沿着边寻找最优解
	if (DEBUG_LOG_TYPE == DEBUG_NOT_INVERTIBLE)
	{
		// 获得顶点坐标的齐次表示
		Vector V1 = {v1->pos.x, v1->pos.y, v1->pos.z, 1};
		Vector V2 = {v2->pos.x, v2->pos.y, v2->pos.z, 1};

		double e1 = dot(V1, Q * V1);
		double e2 = dot(V2, Q * V2);
		double e3 = -1;

		// 沿着边，点可以表示为 v = (1-t) * v1 + t * v2，构造关于 t 的二次函数 a * t^2 + b * t + c
		double a = e1 + e2 - 2 * dot(V1, Q * V2);
		double t = -1;

		// 只有 a > 0 时才有最优解，这里需要 a 足够大，防止溢出
		if (a > ACCURACY)
		{
			double b = 2 * (dot(V1, Q * V2) - e1);
			t = -b / (2 * a);
		}

		// 只有当解在线段中才有效
		if (t >= 0 && t <= 1)
		{
			V = V1 * (1 - t) + V2 * t;
			e3 = dot(V, Q * V);
		}

		// 如果解不在线段中，或者误差是负数，则在两个端点和中点之间选择最优解
		if (e3 < 0)
		{
			Vector V3 = (V1 + V2) / 2;

			// 误差取绝对值，因为奇异情形下误差可能很接近零，结果可能为负数
			e1 = fabs(e1);
			e2 = fabs(e2);
			e3 = fabs(dot(V3, Q * V3));
			V = (e2 < e1) ? ((e3 < e2) ? V3 : V2) : ((e3 < e1) ? V3 : V1);
		}
	}

	// 用最后一位存放误差，因为只有前 3 位表示坐标
	V[3] = fabs(dot(V, Q * V));

	return V;
}

// 曲面简化算法 QEM，不考虑非流形网格：有效对总是边
bool Mesh::Simplify(int times)
{
	// 首先清除误差矩阵列；由于我们会修改顶点位置，因此顶点法向需要清空
	m_nList.clear();
	m_QList.clear();

	// 首先要预处理获得所有有效对。由于不考虑非流形，有效对就是所有半边，什么也不用做
	// 然后对每个顶点计算一个 4x4 误差矩阵 Q，保存一个解向量
	for (int i = 0; i < VertexNum(); ++i)
	{
		vertex3f *v = m_vList[i];
		m_QList.push_back(QEMatrix(v));

		// 对所有顶点标号，方便与 Q 对应
		v->identity = i;
	}

	// 使用红黑树保存半边标记和对应的误差，用一个列表存放有效对的节点指针
	RBT<valid_pair> rbt;
	std::deque<RBT_Node<valid_pair> *> vpList;

	// 循环半边，每个半边就是一对点，计算每个半边上的最小误差
	for (int i = 0; i < m_eList.size(); i++)
	{
		half_edge *e = m_eList[i];

		// 检查是否存在开放边界，如果有就退出
		if (e->twin == nullptr)
			return false;

		// 标记半边，方便与 vpList 中的节点对应
		e->identity = i;

		// 获得最优解，计算此时的误差，将其保存；同时保存插入节点指针
		vpList.push_back(rbt.insert(valid_pair(e, MinVector(e))));
	}

	// 循环 times 次进行简化，剩余面不能少于 4 个
	while (times > 0 && FaceNum() > 4)
	{
		// 从树中得到具有最小误差的顶点合并
		valid_pair vp = rbt.min();
		half_edge *e = vp.e;
		vertex3f *v1 = e->v;
		vertex3f *v2 = Destination(e);

		// 简化操作会删除两个面，这里将两个面上的半边节点移除
		half_edge *e1 = e;
		half_edge *e2 = e->twin;
		for (int i = 0; i < 3; i++)
		{
			rbt.remove(vpList[e1->identity]);
			vpList[e1->identity] = nullptr;
			e1 = e1->next;

			rbt.remove(vpList[e2->identity]);
			vpList[e2->identity] = nullptr;
			e2 = e2->next;
		}

		// 收缩此半边，如果失败，就退出循环；成功则继续
		DEBUG_USER_ASSERT("Simplify", DeleteHalfEdge(e), "Cannot delete half edge.", break)

		// 新顶点直接替换掉对应半边的起点 v1，移动该顶点位置即可
		v1->pos.x = vp.v[0];
		v1->pos.y = vp.v[1];
		v1->pos.z = vp.v[2];

		// 获取以 v1 为起点的所有半边，重新计算顶点误差矩阵和半边误差
		m_QList[v1->identity] = QEMatrix(v1);
		std::deque<half_edge *> eList = AdjacentHalfEdge(v1);
		for (int i = 0; i < eList.size(); ++i)
		{
			half_edge *e = eList[i];

			// 重新计算顶点的误差矩阵
			vertex3f *v2 = Destination(e);
			m_QList[v2->identity] = QEMatrix(v2);

			// 获得最优解
			Vector v = MinVector(e);

			// 注意要先删除所有半边对应 rbt 中的节点，然后重新插入误差
			rbt.remove(vpList[e->identity]);
			rbt.remove(vpList[e->twin->identity]);

			vpList[e->identity] = rbt.insert(valid_pair(e, v));
			vpList[e->twin->identity] = rbt.insert(valid_pair(e->twin, v));
		}
		times--;
	}

	// 再移除一次 QList
	m_QList.clear();

	// 移除所有离散元素
	DeleteIsolated(CGElement::vertex);
	DeleteIsolated(CGElement::half_edge);
	DeleteIsolated(CGElement::face);
	return true;
}

// 设置网格材质
void Mesh::SetMaterial(double M[][3], double alpha)
{
	// 光照系数
	for (int i = 0; i < 3; i++)
		for (int j = 0; j < 3; j++)
			m_material.R[i][j] = M[i][j];

	// 镜面反射系数
	m_material.alpha = alpha;
}

// 网格保存
bool Mesh::SaveMesh(const std::string &filename)
{
	// 根据保存名称修改网格名
	m_name = filename;

	// 首先输出各个属性的描述
	std::fstream fp(filename, std::ios::out);
	fp << "VertexNum " << VertexNum() << std::endl;
	fp << "FaceNum " << FaceNum() << std::endl;
	fp << "HasNormal " << HasVertexNormal() << std::endl;

	// 输出材质信息
	fp << "Material " << m_material.alpha << " ";
	for (int i = 0; i < 3; i++)
		for (int j = 0; j < 3; j++)
			fp << m_material.R[i][j] << " ";
	fp << std::endl;

	// 然后输出顶点信息，并对顶点标号
	for (int i = 0; i < VertexNum(); i++)
	{
		point3f &pos = m_vList[i]->pos;
		m_vList[i]->identity = i;
		fp << "v " << pos.x << " " << pos.y << " " << pos.z << std::endl;
	}
	// 输出面信息
	for (int i = 0; i < FaceNum(); i++)
	{
		fp << "f ";
		std::deque<vertex3f *> vList = VertexOnFace(m_fList[i]);

		// 先输出顶点数，然后输出顶点的索引
		fp << vList.size() << " ";
		for (int j = 0; j < vList.size(); j++)
			fp << vList[j]->identity << " ";

		// 对面编号
		m_fList[i]->identity = i;
		fp << std::endl;
	}
	// 输出邻面信息
	for (int i = 0; i < FaceNum(); i++)
	{
		fp << "a ";

		// 然后我们输出相邻面的索引
		std::deque<face *> fList = AdjacentFace(m_fList[i]);
		fp << fList.size() << " ";
		for (int j = 0; j < fList.size(); j++)
			fp << fList[j]->identity << " ";

		fp << std::endl;
	}

	// 如果有顶点法向就输出顶点法向
	if (HasVertexNormal())
	{
		for (int i = 0; i < VertexNum(); i++)
		{
			vector3f &n = m_nList[i];
			fp << "n " << n.x << " " << n.y << " " << n.z << std::endl;
		}
	}

	// 最后输出终止字符
	fp << "finish" << std::endl;

	return true;
}

// 网格读取
bool Mesh::LoadMesh(const std::string &filename, bool testSame)
{
	// 文件打不开报错
	std::fstream fp(filename, std::ios::in);
	if (!fp.is_open())
		return false;

	// 网格名就是文件名
	m_name = filename;

	std::string detail;
	double value;

	// 首先跳过各个属性的描述
	fp >> detail >> value;
	fp >> detail >> value;
	fp >> detail >> value;

	// 读取材质信息
	fp >> detail >> value;
	double M[3][3];
	for (int i = 0; i < 3; i++)
		for (int j = 0; j < 3; j++)
			fp >> M[i][j];
	SetMaterial(M, value);

	// 暂存顶点，用于重合顶点的情形
	std::deque<vertex3f *> vList;

	// 然后根据第一个字符判断读取内容
	fp >> detail;
	while (detail == "v")
	{
		double x, y, z;
		fp >> x >> y >> z;

		// 禁止重合顶点的检测
		vertex3f *v = CreateVertex({x, y, z}, testSame);
		if (testSame)
			vList.push_back(v);
		fp >> detail;
	}
	// 长度为 6 就足够了
	vertex3f *v[6];
	while (detail == "f")
	{
		int i, n;
		fp >> n;

		// 检查长度
		DEBUG_ASSERT(n <= 6, DEBUG_OUT_OF_RANGE, false)

		// 当开启检测时，使用记录列表中的顶点
		for (int j = 0; j < n; j++)
		{
			fp >> i;
			v[j] = testSame ? vList[i] : m_vList[i];
		}

		// 由于后面会进行衔接操作，因此创建面时取消自动衔接
		CreateFace(v, n, false);
		fp >> detail;
	}
	int i = 0;
	while (detail == "a")
	{
		int size, j;

		// 获得邻面索引信息，然后对这些面进行衔接
		fp >> size;
		for (int k = 0; k < size; k++)
		{
			fp >> j;
			StitchFace(m_fList[i], m_fList[j]);
		}
		fp >> detail;
		i++;
	}
	while (detail == "n")
	{
		double x, y, z;
		fp >> x >> y >> z;
		m_nList.push_back(vector3f(x, y, z));
		fp >> detail;
	}

	return true;
}

// 尝试用网格点创建四面体，返回是否成功
bool CG::tetrahedron(Mesh *mesh)
{
	// 总顶点数至少是 4，并且还没有创建任何面和边
	if (mesh == nullptr || mesh->m_vList.size() < 4 || mesh->m_fList.size() > 0 || mesh->m_eList.size() > 0)
		return false;

	// 先用前 3 个顶点建立一个三角形
	auto it = mesh->m_vList.begin();
	vertex3f *v1 = *it++;
	vertex3f *v2 = *it++;
	vertex3f *v3 = nullptr;
	vertex3f *v4 = nullptr;

	// 第三个顶点不能与前两个共线
	for (; it != mesh->m_vList.end(); ++it)
	{
		// 找到不共线的顶点，方便起见，将它插入到头部
		if (!onLine(v1->pos, v2->pos, (*it)->pos))
		{
			v3 = *it;
			mesh->m_vList.erase(it);
			mesh->m_vList.push_front(v3);
			break;
		}
	}
	// 如果没找到，说明全共线，退出
	if (v3 == nullptr)
		return false;

	// 创建三角形
	face *f = mesh->CreateFace(v1, v2, v3);

	// 寻找第 4 个顶点，不与前三个共面
	it = mesh->m_vList.begin();
	it++, it++, it++;
	for (; it != mesh->m_vList.end(); ++it)
	{
		// 找到不共面的顶点，方便起见，将它插入到头部
		if (!mesh->OnFace((*it)->pos, f))
		{
			v4 = *it;
			mesh->m_vList.erase(it);
			mesh->m_vList.push_front(v4);
			break;
		}
	}
	// 如果没找到，说明全共面，退出
	if (v4 == nullptr)
		return false;

	// 检查 4 顶点是否在面的上方
	bool above = mesh->AboveFace(v4->pos, f);
	if (above)
	{
		// 如果在上方，就要删除原来的面，将顺序反转，以确保法向正确
		mesh->DeleteFace(f);
		mesh->CreateFace(v3, v2, v1);
		mesh->CreateFace(v1, v4, v3);
		mesh->CreateFace(v3, v4, v2);
		mesh->CreateFace(v1, v2, v4);
	}
	else
	{
		mesh->CreateFace(v1, v3, v4);
		mesh->CreateFace(v3, v2, v4);
		mesh->CreateFace(v1, v4, v2);
	}

	// 最后要移除所有孤立面
	mesh->DeleteIsolated(CGElement::face);

	return true;
}

// 取网格点的凸包
bool CG::convexHull(Mesh *mesh)
{
	// 先尝试创建一个四面体
	if (!tetrahedron(mesh))
		return false;

	// 增量算法，对每个四面体以外的顶点，计算被这个顶点照射到的面
	// 将这些面删除，然后把面的边与顶点建立新的面
	auto vit = mesh->m_vList.begin();
	vit++, vit++, vit++, vit++;
	for (; vit != mesh->m_vList.end(); ++vit)
	{
		// 标记被照射的面，只需要保存一个面的指针，因为所有标记面都相邻
		std::deque<face *> fList;
		for (auto fit = mesh->m_fList.begin(); fit != mesh->m_fList.end(); ++fit)
		{
			// 在表面上面，说明被照射到，将其记录下来
			if (mesh->AboveFace((*vit)->pos, *fit))
			{
				(*fit)->identity = 1;
				fList.push_back(*fit);
			}
			else
				(*fit)->identity = 0;
		}

		// 检查被照射的面，标记相邻的未照射面
		std::deque<half_edge *> eList;
		for (auto fit = fList.begin(); fit != fList.end(); fit++)
		{
			// 检查面上的半边，记录边界边，同时将边界边上的顶点保存的半边修改为对边，防止在删除面时孤立
			std::deque<half_edge *> eList0 = mesh->HalfEdgeOnFace(*fit);
			for (int i = 0; i < eList0.size(); i++)
			{
				half_edge *e = eList0[i];
				if (e->twin->incFace->identity != 1)
				{
					mesh->Destination(e)->e = e->twin;
					eList.push_back(e->twin);
				}
			}
		}

		// 删除光照面。直接删除不会有问题，如果出现顶点指向的边为空的情况
		// 有两种可能：1. 离散点，则在最后会被删除；2. 处于被光照和未被光照的边界上，这种情况在前面进行修正
		for (auto it = fList.begin(); it != fList.end(); ++it)
			mesh->DeleteFace(*it);

		// 然后去除离散的面
		mesh->DeleteIsolated(CGElement::face);

		// 遍历所有反向边，按照连接顺序，每次要创建具有公共边的新面，防止出现非流形网格
		// 这样创建的面都会相邻，不会出现恰好有一个公共顶点的两个面，而周围都是空的情况
		int count = eList.size();
		auto eit = eList.begin();
		while (count > 0)
		{
			mesh->CreateFace(mesh->Destination((*eit)), (*eit)->v, *vit);
			count--;

			// 找到目标点对应的边
			for (auto it = eList.begin(); it != eList.end(); ++it)
			{
				if ((*it)->v == mesh->Destination((*eit)))
				{
					eit = it;
					break;
				}
			}
		}
	}
	// 最后移除所有离散点
	mesh->DeleteIsolated(CGElement::vertex);
	return true;
}

// 检查点是否在凸包内部
bool CG::isInHull(point3f p, std::deque<point3f> &hull)
{
	// 利用二分角检测
	int n = hull.size();
	int a = 1, b = n - 1;
	vector3f p0 = p - hull[0];

	// 检查点的顺序是逆时针还是顺时针
	vector3f pa = hull[a] - hull[0];
	vector3f pb = hull[b] - hull[0];

	// 叉积规范化作为法向
	vector3f normal = normalize(cross(pa, pb));

	// 循环直到两边界点相邻
	while (abs(a - b) > 1)
	{
		// 获得等分点与初始顶点的向量
		int mid = (a + b) / 2;
		vector3f m0 = hull[mid] - hull[0];
		double sign = det(p0, m0, normal);

		// 判断如果恰好在连线上，就计算距离，如果距离满足直线关系，就说明在内部
		if (fabs(sign) < ACCURACY)
		{
			if (fabs(distance(hull[0], p) + distance(hull[mid], p) - distance(hull[0], hull[mid])) < ACCURACY)
				return true;
			else
				return false;
		}

		// 检查点在等分线的哪一侧，同时划分多边形
		if (sign > 0)
			b = mid;
		else
			a = mid;
	}

	// 剩下的一定是三角形，只要判断是否在三角形内部
	vector3f p0pa = cross(hull[0] - p, hull[a] - p);
	vector3f papb = cross(hull[a] - p, hull[b] - p);
	vector3f pbp0 = cross(hull[b] - p, hull[0] - p);

	// 如果向量方向不同，说明在外部
	int s1 = SIGN(dot(p0pa, papb));
	int s2 = SIGN(dot(papb, pbp0));
	int s3 = SIGN(dot(pbp0, p0pa));
	if (s1 == s2 && s2 == s3)
		return true;

	return false;
}

// 计算并返回平面多边形的凸包
std::deque<point3f> CG::convexHull(std::deque<point3f> &polygon)
{
	// 按照顺序检查相邻两个边向量的叉积的符号，如果出现符号反转，说明两个边的公共顶点“凹陷“，将异号的顶点移除后得到凸包
	std::deque<point3f> hull = polygon;

	// 三角形一定是凸包，作为退出条件
	while (hull.size() > 3)
	{
		int n = hull.size();

		std::deque<vector3f> pro(n);
		for (int i = 0; i < n; i++)
		{
			// 获得前后两个顶点的索引
			int indp1 = (i - 1 + n) % n;
			int indp2 = (i + 1) % n;
			vector3f p1 = hull[indp1] - hull[i];
			vector3f p2 = hull[i] - hull[indp2];

			// 保存叉积结果
			pro[i] = cross(p1, p2);
		}

		// 检查这些叉积的方向，同时统计符号
		int pos = 1;
		std::deque<int> index(n);
		for (int i = 1; i < n; i++)
		{
			// 全与第一个叉积做点积获得符号
			index[i] = SIGN(dot(pro[i], pro[0]));
			if (index[i] > 0)
				pos++;
		}

		// 如果全部同号，说明已经形成凸包
		if (pos == n)
			break;

		// 删除对应少数的符号
		int sign = (pos > (n - pos)) ? -1 : 1;
		auto it = hull.begin();
		for (int i = 0; i < n; i++)
		{
			if (index[i] == sign)
				it = hull.erase(it);
			else
				it++;
		}
	}

	return hull;
}

// 球面：中心位置、三个半长轴，分割份数
bool CG::sphere(Mesh *mesh, point3f center, vector3f axis, double a, double b, double c, int div)
{
	// 首先检查参数有效性，分割份数至少应该是 2
	if (mesh == nullptr || div < 1)
		return false;

	// 把 a,b,c 都换成正的
	a = fabs(a);
	b = fabs(b);
	c = fabs(c);

	// 使用极坐标取点，获得纵向步长，也是横向步长
	const double PI = 3.1415926;
	double dphi = PI / div;
	double phi = dphi;
	double theta = 0;

	// 获得从 z 轴转动到 axis 的四元数
	axis = normalize(axis);
	quaternion q = rotate({0, 0, 1}, axis);

	// 不创建极点，i 是纵向，因此纵向点的数量少一些，j 是横向，横向分割份数是纵向的 2 倍
	for (int i = 0; i < div - 1; i++)
	{
		for (int j = 0; j < div * 2; j++)
		{
			// 先记录一个纯四元数，将其转动
			quaternion p({a * sin(phi) * cos(theta), b * sin(phi) * sin(theta), c * cos(phi)});
			p = q * p / q;
			mesh->CreateVertex({p.v.x + center.x, p.v.y + center.y, p.v.z + center.z});
			theta += dphi;
		}
		theta = 0;
		phi += dphi;
	}

	// 现在创建这些点构成的面，每一行和上一行构成面
	for (int i = 1; i < div - 1; i++)
	{
		for (int j = 0; j < div * 2; j++)
		{
			vertex3f *v0 = mesh->m_vList[i * (div * 2) + j];
			vertex3f *v1 = mesh->m_vList[i * (div * 2) + (j + 1) % (div * 2)];
			vertex3f *v2 = mesh->m_vList[(i - 1) * (div * 2) + (j + 1) % (div * 2)];
			vertex3f *v3 = mesh->m_vList[(i - 1) * (div * 2) + j];
			mesh->CreateFace(v0, v1, v2, v3);
		}
	}

	// 最后创建极点，连接极点和相邻的点
	point3f up = center + axis * c;
	mesh->CreateVertex(up);
	for (int j = 0; j < div * 2; j++)
	{
		vertex3f *v0 = mesh->m_vList[j];
		vertex3f *v1 = mesh->m_vList[(j + 1) % (div * 2)];
		vertex3f *v2 = mesh->m_vList.back();
		mesh->CreateFace(v0, v1, v2);
	}

	point3f dp = center + (-axis) * c;
	mesh->CreateVertex(dp);
	for (int j = 0; j < div * 2; j++)
	{
		vertex3f *v0 = mesh->m_vList[(div - 2) * (div * 2) + j];
		vertex3f *v1 = mesh->m_vList[(div - 2) * (div * 2) + (j + 1) % (div * 2)];
		vertex3f *v2 = mesh->m_vList.back();
		mesh->CreateFace(v1, v0, v2);
	}
	return true;
}

// 柱面：中心位置，中轴，半径，高度，分割份数
bool CG::cylinder(Mesh *mesh, point3f center, vector3f axis, double r, double h, int divr, int divh)
{
	// 首先检查参数有效性，分割份数至少应该是 2
	if (mesh == nullptr || divr < 1 || divh < 1)
		return false;

	// 把 r,h 都换成正的
	r = fabs(r);
	h = fabs(h);

	// 使用极坐标取点，获得纵向步长，和横向步长
	const double PI = 3.1415926;
	double dtheta = 2 * PI / divr;
	double theta = 0;
	double dh = h / divh;

	// 获得从 z 轴转动到 axis 的四元数
	quaternion q = rotate({0, 0, 1}, axis);

	double nz = -h / 2;
	// 从下向上依次创建顶点
	for (int i = 0; i < divh + 1; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			// 先记录一个纯四元数，将其转动
			quaternion p({r * cos(theta), r * sin(theta), nz});
			p = q * p / q;
			mesh->CreateVertex({p.v.x + center.x, p.v.y + center.y, p.v.z + center.z});
			theta += dtheta;
		}
		theta = 0;
		nz += dh;
	}

	// 现在创建这些点构成的面，每一行和上一行构成面
	for (int i = 0; i < divh; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			vertex3f *v0 = mesh->m_vList[i * divr + j];
			vertex3f *v1 = mesh->m_vList[i * divr + (j + 1) % divr];
			vertex3f *v2 = mesh->m_vList[(i + 1) * divr + (j + 1) % divr];
			vertex3f *v3 = mesh->m_vList[(i + 1) * divr + j];
			mesh->CreateFace(v0, v1, v2, v3);
		}
	}

	// 创建上下两个面的中点
	quaternion p1({0, 0, -h / 2});
	quaternion p2({0, 0, h / 2});
	p1 = q * p1 / q;
	p2 = q * p2 / q;
	vertex3f *c1 = mesh->CreateVertex({p1.v.x + center.x, p1.v.y + center.y, p1.v.z + center.z});
	vertex3f *c2 = mesh->CreateVertex({p2.v.x + center.x, p2.v.y + center.y, p2.v.z + center.z});

	// 建立上下两个面
	for (int i = 0; i < divr; i++)
	{
		vertex3f *v0 = mesh->m_vList[i];
		vertex3f *v1 = mesh->m_vList[(i + 1) % divr];
		mesh->CreateFace(c1, v1, v0);

		vertex3f *v2 = mesh->m_vList[divh * divr + i];
		vertex3f *v3 = mesh->m_vList[divh * divr + (i + 1) % divr];
		mesh->CreateFace(v2, v3, c2);
	}

	return true;
}

// 锥面：中心位置，中轴，半角，高度，分割份数
bool CG::cone(Mesh *mesh, point3f origin, vector3f axis, double theta, double h, int divr, int divh)
{
	// 首先检查参数有效性，分割份数至少应该是 2
	if (mesh == nullptr || theta <= 0 || divr < 1 || divh < 1)
		return false;

	// 把 h 换成正的，theta 转换为弧度
	const double PI = 3.1415926;
	h = fabs(h);
	theta = theta / 180 * PI;

	// 使用极坐标取点，获得纵向步长，和横向步长
	double dphi = 2 * PI / divr;
	double phi = 0;
	double dh = h / divh;

	// 获得从 z 轴转动到 axis 的四元数
	quaternion q = rotate({0, 0, 1}, axis);

	// 先跳过极点
	double nz = dh;
	// 从下向上依次创建顶点
	for (int i = 0; i < divh; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			// 先记录一个纯四元数，将其转动
			double r = nz * tan(theta);
			quaternion p({r * cos(phi), r * sin(phi), nz});
			p = q * p / q;
			mesh->CreateVertex({p.v.x + origin.x, p.v.y + origin.y, p.v.z + origin.z});
			phi += dphi;
		}
		phi = 0;
		nz += dh;
	}

	// 现在创建这些点构成的面，每一行和上一行构成面
	for (int i = 0; i < divh - 1; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			vertex3f *v0 = mesh->m_vList[i * divr + j];
			vertex3f *v1 = mesh->m_vList[i * divr + (j + 1) % divr];
			vertex3f *v2 = mesh->m_vList[(i + 1) * divr + (j + 1) % divr];
			vertex3f *v3 = mesh->m_vList[(i + 1) * divr + j];
			mesh->CreateFace(v0, v1, v2, v3);
		}
	}

	// 最后创建锥顶，连接第一行
	mesh->CreateVertex(origin);
	for (int j = 0; j < divr; j++)
	{
		vertex3f *v0 = mesh->m_vList[j];
		vertex3f *v1 = mesh->m_vList[(j + 1) % divr];
		vertex3f *v2 = mesh->m_vList.back();
		mesh->CreateFace(v2, v1, v0);
	}

	// 创建下面的中点
	quaternion p({0, 0, h});
	p = q * p / q;
	mesh->CreateVertex({p.v.x + origin.x, p.v.y + origin.y, p.v.z + origin.z});

	// 建立下面
	for (int j = 0; j < divr; j++)
	{
		vertex3f *v0 = mesh->m_vList[(divh - 1) * divr + j];
		vertex3f *v1 = mesh->m_vList[(divh - 1) * divr + (j + 1) % divr];
		vertex3f *v2 = mesh->m_vList.back();
		mesh->CreateFace(v0, v1, v2);
	}

	return true;
}

// 圆环面：中心位置，垂直轴，内径，外径，分割份数
bool CG::ring(Mesh *mesh, point3f center, vector3f axis, double r, double R, int divR, int divr)
{
	// 首先检查参数有效性，分割份数至少应该是 2
	if (mesh == nullptr || divr < 1 || divR < 1 || fabs(R) <= fabs(r))
		return false;

	// 把 R,r 换成正的
	const double PI = 3.1415926;
	R = fabs(R);
	r = fabs(r);

	// 使用极坐标取点，大圆步长和小圆步长
	double dRa = 2 * PI / divR;
	double dra = 2 * PI / divr;
	double Ra = 0;
	double ra = 0;
	double mR = (R + r) / 2;
	double mr = (R - r) / 2;

	// 获得从 z 轴转动到 axis 的四元数
	quaternion q = rotate({0, 0, 1}, axis);

	// 转动创建顶点
	for (int i = 0; i < divR; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			// 先记录一个纯四元数，将其转动
			quaternion p({(mR + mr * sin(ra)) * cos(Ra), (mR + mr * sin(ra)) * sin(Ra), mr * cos(ra)});
			p = q * p / q;
			mesh->CreateVertex({p.v.x + center.x, p.v.y + center.y, p.v.z + center.z});
			ra += dra;
		}
		ra = 0;
		Ra += dRa;
	}

	// 现在创建这些点构成的面，每一行和上一行构成面
	for (int i = 0; i < divR; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			vertex3f *v0 = mesh->m_vList[i * divr + j];
			vertex3f *v1 = mesh->m_vList[i * divr + (j + 1) % divr];
			vertex3f *v2 = mesh->m_vList[((i + 1) % divR) * divr + (j + 1) % divr];
			vertex3f *v3 = mesh->m_vList[((i + 1) % divR) * divr + j];
			mesh->CreateFace(v0, v1, v2, v3);
		}
	}
	return true;
}